##
## Distance to powerplant
##

options(rasterTmpDir='C:/temp/')
require(Rcpp)

setwd(wkdir)
output_version <- Sys.Date()

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01',crs=4326)

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates() 

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones = st_drop_geometry(subset(ntl.zones,select=c(OBJECTID, GID_0)))
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,ntl.zones)

nCpu = 6

##
## Powerplants AICD
## 

##dist2airafr     double  %9.0g                 Distance from grid pt to any airport in Africa
##dist2airlcd     double  %9.0g                 Distance from grid pt to any airport in the 4 lake chad countries
##dist2airctr     float   %9.0g                 Distance from grid pt to any airport within the country

ntl.zones.pts$dist2ppctr<-NA

# CMR
pp_cmr <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/CMR_PowerPlants.shp'))
pp_cmr <- subset(pp_cmr,STATUS=="OPR")
st_crs(pp_cmr) <- st_crs(ntl.zones.pts)
ntl.zones.pts$dist2ppctr[ntl.zones.pts$GID_0=="CMR"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="CMR",], 
                                                                            pp_cmr, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

# NER
pp_ner <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/NER_PowerPlants.shp'))
pp_ner <- subset(pp_ner,STATUS=="OPR")
st_crs(pp_ner) <- st_crs(ntl.zones.pts)
ntl.zones.pts$dist2ppctr[ntl.zones.pts$GID_0=="NER"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NER",], 
                                                                            pp_ner, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

# NGA
pp_nga <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/NGA_PowerPlants.shp'))
pp_nga <- subset(pp_nga,STATUS=="OPR")
st_crs(pp_nga) <- st_crs(ntl.zones.pts)
ntl.zones.pts$dist2ppctr[ntl.zones.pts$GID_0=="NGA"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="NGA",], 
                                                                            pp_nga, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

# TCD
pp_tcd <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/TCD_PowerPlants.shp'))
pp_tcd <- subset(pp_tcd,STATUS=="OPR")
st_crs(pp_tcd) <- st_crs(ntl.zones.pts)
ntl.zones.pts$dist2ppctr[ntl.zones.pts$GID_0=="TCD"] <- unlist(nngeo::st_nn(ntl.zones.pts[ntl.zones.pts$GID_0=="TCD",], 
                                                                            pp_nga, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000
## 4 countries 
pp_lcd <- rbind(pp_cmr,pp_ner,pp_nga,pp_tcd)
st_crs(pp_lcd) <- st_crs(ntl.zones.pts)
ntl.zones.pts$dist2pplcd<-NA
ntl.zones.pts$dist2pplcd <- unlist(nngeo::st_nn(ntl.zones.pts, pp_lcd, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

## 4 countries + 10 neighbours
pp_bfa <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/BFA_PowerPlants.shp'),crs=4326)
pp_caf <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/CAF_PowerPlants.shp'),crs=4326)
pp_cog <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/COG_PowerPlants.shp'),crs=4326)
pp_gab <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/GAB_PowerPlants.shp'),crs=4326)
pp_gnq <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/GNQ_PowerPlants.shp'),crs=4326)
pp_mli <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/MLI_PowerPlants.shp'),crs=4326)
pp_sdn <- st_read(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/powerplants_aicd/SDN_PowerPlants.shp'),crs=4326)

pp_afr <- rbind(pp_bfa,pp_caf,pp_cmr,pp_cog,pp_gab,pp_gnq,pp_mli,pp_ner,pp_nga,pp_sdn,pp_tcd)
pp_afr <- subset(pp_afr,STATUS=="OPR")

ntl.zones.pts$dist2ppafr<-NA
ntl.zones.pts$dist2ppafr <- unlist(nngeo::st_nn(ntl.zones.pts, pp_afr, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000

##
## WRI 
## Global Energy Observatory, Google, KTH Royal Institute of Technology in Stockholm, Enipedia, World Resources Institute. 
## 2018. Global Power Plant Database. Published on Resource Watch and Google Earth Engine; 
## http://resourcewatch.org/ https://earthengine.google.com/
##
ppwri <- fread(paste0(wkdir,'/local/gis_data/019_utilitiesCommunication/PowerPlants_WRI/global_power_plant_database.csv'))
names(ppwri)
unique(ppwri$country)

ppwri_sf <- st_as_sf(ppwri, coords = c("longitude", "latitude"), crs = 4326)

## all africa
ntl.zones.pts$dist2ppwriafr<-NA
ntl.zones.pts$dist2ppwriafr <- unlist(nngeo::st_nn(ntl.zones.pts, ppwri_sf, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000
summary(ntl.zones.pts)

## region
ppwri_sf <- subset(ppwri_sf,country %in% c('CMR','NER','NGA','TCD'))

ntl.zones.pts$dist2ppwrilcd<-NA
ntl.zones.pts$dist2ppwrilcd <- unlist(nngeo::st_nn(ntl.zones.pts, ppwri_sf, k = 1, returnDist = T, parallel = nCpu)$dist) / 1000
summary(ntl.zones.pts)

save.dta13(st_drop_geometry(ntl.zones.pts), paste0(wkdir,"/proc_data/DIST2PP_all",output_version,".dta"))